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Abstract. We investigate autogenous fragmentation of dry granular material in rotating cylinders using 
two-dimensional molecular dynamics. By evaluation of spatial force distributions achieved numerically for 
various rotation velocities we argue that comminution occurs mainly due to the existence of force chains. A 
statistical analysis of theses force chains explains the spatial distribution of comminution efficiency in ball 
mil ls as measured experimentally by Rothkegel JTJ and Rolf [Bj. For animated sequ ences of our simulations 
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see http:/ /summa. physik. hu-berlin.de/~kies/Research/RotatingCylinder/rotatingcylinder. html 



PACS. 81.05.Rm Porous materials; granular materials - 83.70.-n Granular systems 



1 Introduction 



In the past years molecular dynamics simulations of gran- 
ular systems - such as sand, fertilizer, grain and others 
- have been developed towards a reliable tool for an in- 
vestigation of granular systems. Unlike a few years ago, 
when physicists started to investigate granular assemblies 
consisting of just a few hundred spherical granular parti- 
cles in two dimensions, e.g. ^rfPflR FBl , and less than 
100 particles in three dimensions, e.g. |10|, today, owing 
to the rapid evolution of computer facilities, we are able 
to simulate complex systems of many thousands of par- 
ticles in two and three dimensions even accounting for 
non-sphcrical grain features. 

Due to this progress molecular dynamics of granular 
systems by now can be applied to the simulation of techno- 
logically important processes striving for a better under- 
standing of relevant details of the process and, henceforth, 
for an optimization of industrial technologies (e.g. [Tl| , p^ , [T3| ,| 
Molecular simulation techniques, therefore, might develop 
to a prospective engineering tool partly replacing costly 
laboratory experiments by "computer experiments" . 

In the present paper we apply the method of molecular 
dynamics to the simulation of autogenous dry comminu- 
tion in a tumbling ball mill |jq] . To this end we devised 
a novel algorithm which accounts for the fragmentation 
of particles. First we will demonstrate that experimen- 
tally known results, in particular those by Rothkegel and 
Rolf jjJH , can be reproduced up to a good degree of ac- 
curacy. We will show that the efficiency of comminution 
in ball mills is mainly determined by the presence of force 
chains. Our result will hint at how to improve the effi- 
ciency of a milling machinery widely applied in industry. 



2 The Simulation Model 

2.1 Molecular Dynamics 

For the simulation we assumed pairwise interaction forces 
and applied the model by Cundall and Strack Q and Haff 
and Werner ||: A pair of two-dimensional spherical parti- 
cles i and j interacts through the force 



(1) 



if their distance is smaller than the sum of their radii 



Ri + Rj 



£ii > 0. The unit vectors in normal 



and tangential direction are defined as 



n = r ij /|r ij | 
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TO 3 |lffij l7j^-r, — Tj, The forces in normal and tangential 
direction are respectively given by 
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Here Ri, rrij, rj, and fli are respectively the radius, 
the mass, the position, the velocity and the rotation ve- 
locity of the ith particle. Y = 8 ■ 10 6 gs~ 2 is the elastic 
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constant and 7^ = 800s" 1 , jt = 3000s _1 are the param- 
eters of damping in normal and tangential direction. The 
Coulomb friction constant was assumed to be /i = 0.5. 
These empirical values which we used throughout all sim- 
ulations have proven to render realistic behaviour for a 
typical granular material. 

The normal force (Q) is composed of an elastic repulsive 
part and a dissipative part which acts against the direc- 
tion of motion. For a collision of two-dimensional spheres, 
i.e. disks, the Hertz contact law ~ £J 2 p^ , p0| reduces 

to |3 

where E r is the reduced elastic module, i.e. a material 
constant. Equation (||) provides a relation £ = ^ (F w ), 
however, in the molecular dynamics simulations we need 
the inverse 

f n _ p N (£). To calculate the forces due to 
Eq. (||) one would have to invert (||) numerically for each 
particle contact in each time step, i.e., to solve a tran- 
scendental equation. Alternatively one could tabulate the 
function F N but due to large force gradients, in par- 
ticular in the instant of contact, we have to simulate with 
double precision accuracy. Therefore, the table would be 
extremely large. Apart from the small logarithmic term 
— F?? In F^ in Eq. (Q) one finds that the force is propor- 
tional to the compression £y . Therefore, in our simulations 
we used the linear law (Q) which is a good approximation 
in the force interval of interest. 

Equation (||) describes the relative velocity of the par- 
ticle surfaces at the point of contact which results from 
both the tangential part of the relative particle velocity 
and the velocity of particle spin. 

The Coulomb friction law is taken into account in 
Eq. (||): If the tangential force of two colliding particles 
exceeds fj, times the normal force the particles slide upon 
each other feeling constant friction. In this way the Coulomb 
law formulates a maximum transferable shear force. 

As numerical integration scheme we used the Gear 
predictor-corrector method of sixth order (e.g. |2^| ). 

Let us remark that there exist various models with 
diverse descriptions of the interaction forces, e.g. In 
three-dimensional calculations it is crucial that the dissi- 
pative term in (|J) is proportional to [||) (in viscoelas- 
tic approximation) . In three dimensions the used term ~ £ 
in Eq. (@) yields wrong results |||||. 

Another class of molecular dynamics are event-driven 
simulations, e.g. p6f which require much less numerical ef- 
fort than the method described above (force-driven Molec- 
ular Dynamics). In event-driven simulations one does not 
integrate the equations of motion explicitly but rather ex- 
presses the loss of mechanical energy in normal and tan- 
gential direction due to a collision via coefficients of resti- 
tution e N and e T both of which can be determined from 
the material properties of the particles |24|,|27]]. Besides 
some conceptional difficulties which can be overcome by 
applying numerical tricks (2^] these algorithms are only 
applicable in case the duration of collisions is negligible in 



comparison with the mean free flight time, i.e. the time 
a particle moves without interaction. In this limit colli- 
sions can be assumed as instantaneous events. Moreover, 
this implies that three-particle interactions are very rare 
events. However, in the system under study, namely the 
ball mill, this premise docs not hold because the majority 
of particles permanently is in close contact with neigh- 
bouring particles. 

2.2 Fragmentation Probability of Grains 

It can be observed in experiments that different particles 
of identical size and material vary with respect to their 
tensile strengths. That means when compressing particles 
of a sample by applying a fixed force only a fraction of the 
sample actually will break. For a sufficiently large sample 
this fraction can be rephrased in terms of a fracture prob- 
ability. 

The reason for this different behavior can be explained 
by the existence of flaws which provide sites for stress 
concentration and the initiation of a crack which, subse- 
quently, propagates through the material. These flaws are 
distributed throughout the whole volume, however, sev- 
eral investigations have revealed that surface flaws acti- 
vated by high tensile stress play the dominant role for the 
initiation of a crack [^9). Hence, the fracture probability 
is intimately related to the statistical distribution of sur- 
face flaws and the resulting stress distribution over the 
particle's surface. 

Starting from some plausible assumptions concerning 
the flaw distribution and the behaviour of the material the 
fracture probability P of a single particle can be derived 
employing statistical reasoning. A particle of radius i?, 
subjected to a force which stores the specific elastic energy 
W m = W/m, will break with probability |3C| 

P{R, W m ) ~ 1 ~ exp (~cR 2 W^) . (9) 

Here, c and z are material constants. The original deriva- 
tion (considering rods of a brittle material) goes back to 
Weibull (3^] . This simple law has proven in practice by fit- 
ting various experimental data. The constants c and z can 
be extracted from measurements when plotting [lnln(l — 
P) _1 ] vs. [lnW m ]. Typical values of z are found to lie in 
the range 1.5, . . . , 2.5. 

Notice that in agreement with experimental results the 
term R 2 in the exponent (for fixed W m ) in general predicts 
an anticorrclation between particle size and resistance to 
breakage. Clearly, this fact is rooted in the diminished 
probability to find a sufficiently large flaw on a smaller 
surface. 

For spherical particles (or, more generally, for particles 
having convex surface), assuming linear elastic behaviour, 
Hertz-theory Jl^,0 can be employed to derive a relation 
between the elastic energy W stored in a sphere and the 
exerted repulsive contact force F 

1 

2 2 / D 2 F 5 \ 3 

w -¥ F -l{ £ f) ■ (10 > 
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where h quantifies the particle's deformation and D is 
another material constant (expressed by the Poisson ratio 
and constant of elasticity) |l^,^0). It should be kept in 
mind that formula (10) remains valid in the linear elastic 
regime only. 

Applied to our simulation of ball mills we make use 
of these laws in the following way: Consider the situation 
sketched in Fig. |]: A particle placed in the centre is com- 
pressed by four impacting particles. 




Fig. 1. A typical situation of one particle being compressed 
by four impacting particles. The forces are considered as inde- 
pendent stress sources. Therefore, in simulations we evaluated 
only the maximum of all exerted forces to compute the fracture 
probability P(R,W m ). 



The highest compressive stress occurs around the con- 
tact point; a fact well known to engineers |3^ | and also 
found in numerical simulations |Bjjj34j . Hence, the fracture 
probability of the centre particle is linked to the probabil- 
ity of finding a flaw of sufficient size in the very vicinity 
of any of the four contact regions. Thus, the action of dif- 
ferent impacting particles can be regarded as independent 
stress sources. That means we consider the elastic energy 
stored by each impacting particle separately, calculate the 
related fracture probability according to formulae (^|) and 
( |l(]| ) as a function of the particle's size and the maximum 
of all acting normal forces F N 



P{R,F N ) ~ 1-exp 



-c7?( 2 ~ z ^) (F 



(11) 



and actually decide the question of breakage by subse- 
quently drawing a random number for each of the contact 
sites. 

The constant c is related to the empirical constant c in 
the precise connection is elaborated by inserting 
|) into (^|) and expressing m by the particle radius 

m oc Rr 

Up to this point we discussed the fracture probabil- 
ity for three-dimensional spheres. In our simulations we 
modelled two-dimensional spheres, i.e. circular discs of 
unique thickness L. The fracture probability correspond- 
ing to Eq. (Ill]) can be obtained in analogy to (0) yielding 




Note that the surface of a sphere oc R 2 has to be replaced 
by the surface of a cylinder oc RL (excluding top and 
bottom wall). Correspondingly Eq. ( |l0|) turns into h = 
YF N and, hence, the elastically stored energy for spheres 
has to be replaced by 



W = \hF N = ^ ~ (F N ) 2 



(13) 



Taking into account that the mass of a cylinder is propor- 
tional oc R 2 , one obtains 



P{R, F N ) - 1 - exp ( -cR 



(F N ) 2 



R? 



(14) 



In section || we will motivate the choice z = 2. With this 
specification one finally arrives at the following expression 
for the fracture probability 



P(R, F N )~1- cxp (-ciT 3 (F N ) 



(15) 



3 Fragment size distribution 

Once a flaw has been activated it forms an initial crack 
which rapidly propagates through the particle. Typical 
crack velocities in the range of 1500m/s have been mea- 
sured j?9| . The formation of a variety of differently sized 
fragments occurs by virtue of branching cascades. Branch- 
ing is governed by a balance between the energy release 
rate G and the so-called crack resistivity B [p9| the lat- 
ter accounting for the creation of new surfaces. Because of 
the extremely short duration of the rupture process one 
can safely neglect external energetic contributions to this 
energy balance. The energy release rate G depends on the 
crack velocity and the crack length whereas the crack resis- 
tivity B depends on the crack velocity and on the number 
of branches. Due to a maximal crack velocity, i.e. the speed 
of sound, the balance between both terms requires the for- 
mation of new cracks. From these considerations it can be 
shown J2t| that the fragment size distribution (cumulative 
mass distribution) Q is a function of the product R f W m , 
where Rf denotes the fragment size, i.e. Q = Q(RfW, 

Experimental evidence |3(J for a scaling law Q 
in a gree ment with the well-known empirical Schuhmann 
law [B5l 

k 



R? f ia 



(16) 



which itself can be regarded as an approximation of the 
Rosin-Rammler law [M 



1 — exp 



k 



(17) 



P(R, W m )~l- exp (-cRLW* ) 



(12) 



The variable k has the dimension of a length and, thus, 
can only be identified with the size of the original particle. 

A rigorous derivation of Eq. (17) - compatible with an 
exponent = 1 — starting from rather mild assumption 
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and applying Poisson statistics was performed by Gilvary 
(see also @). 

It has to be mentioned that the derived laws Eqs. ( |l6|Jr7| ) 
describe the distribution of fragments only below a certain 
size (endoclastic vs. exoclastic distribution p7j). For the 
largest fragments an equivalently well-defined quantita- 
tive statement does not exist. The deviation between theo- 
retical predictions and experimental observations typically 
occurs for grains which together collect about 75-90% of 
the total mass ^]. Moreover, some diverging theoretical 
mean values can be understood from the fact that the 
assumption of equidistributed flaws on all length scales 
neglects the effect of flaw depletion applying to tiny frag- 
ments. Flaw depletion thus sets a lower bound to the 
range of applicability. Insofar the size of fragments which 
is satisfactorily described by the distributions Eq. ( |l6| ) and 
Eq. (|j]) is restricted to the intermediate range. 

When fitting both the more refined Rosin-Rammler 
law Eq. ([It]) and the simplified Schumann law Eq. ( |l6| ) to 
diverse experimental data it was found that in some situ- 
ations the Rosin-Rammler law proved superior [p9| (mill 
products P8[l) whereas in other cases (specimen fractured 
in gelatin |38|) the Schuhmann law was even more ade- 
quate. 

Finally, we mention the population balance model p9| 
which assumes that the fragment size distribution can oe 
normalized to the initial particle size, i.e. Q — Q(Rj/R). 
This assumption is based on empirical observations. More- 
over, its central assertion can also be derived when start- 
ing from a Weibull exponent z = 2 and applying statistical 
reasoning similar to the one sketched above. 

In our simulations we employed the Rosin-Rammler 
law Eq. (jl|) with k = R and an exponent (3=1. Since 
experimental data have shown the mean dimension of the 
largest fragment to be roughly 75% of the original particle 
size [[38], i.e. RJ 1 ^ = the normalized fragment size 
distribution simply reads 



responds to a situation of extreme compression an unre- 
alistically strong repulsive normal force together with 
an exceptional deformation energy would be the conse- 
quence. We bypassed this problem by modifying the in- 
teraction between the two fragments in the following way: 
Up to the moment of first complete separation of two such 
fragments i and j the normal force ^ at the fc + l-rst time 
step is replaced by 



rpN _ 



with 



y#,(*)-m: 



y 7*(*t -, 3 

if (r ? ±j 

otherwise 



ij) • n 



• n< 



(19) 



C#) - |r,(fc - 1) - rj(k - 1)| - \n(k) - rj(k)\ . (20) 

This means that the two fragments, on the one hand, re- 
sist to further compression like standard particles and, on 
the other hand, experience a small repelling force e which 
drives them apart. The repelling force is chosen sufficiently 
small in order to suppress unphysical energy gain of the 
system; this means that the unavoidable energy, which due 
to this force is pumped into the system anyhow, must be 
negligible in comparison with the mean kinetic energy of 
the grains. After the moment of first complete separation 
of the two fragments they also interact in the standard 
way, as dictated by Eq. (|l|). Of course, the interaction of 
two such fragments with all other particles is never modi- 
fied. By trimming of e we could not only suppress the elas- 
tic energy gain but also evidenced that the period requir- 
ing a modified interaction force was comparatively short. 
The procedure of fragmentation is sketched in Fig. ||. We 
want to mention that multiple fragmentation, i.e. further 
fragmentation of the fragments is included in the model 
and is frequently observed in simulations. 



1 - exp ( 



1 — exp — 



(18) 



3.1 Molecular Dynamics Modelling of Fragmentation 

In the present study we exclusively utilize spherical par- 
ticles (disks), which means that also the fragments of a 
disrupted particle must be spherical. The single fracture 
statistics, as discussed in section [|, is put into practice by 
dismembering a cracking disk into two fragments. Whereas 
the size of one fragment is chosen in accordance with the 
Rosin-Rammler law the size of the second one is deter- 
mined by mass conservation, i.e. for two-dimensional par- 
ticles by area conservation. 

In the majority of situations a breaking particle is 
closely surrounded by neighbouring ones. Hence, in gen- 
eral, it will be impossible to place the two fragments avoid- 
ing an overlap. However, since a considerable overlap cor- 




Fig. 2. Sketch of the particle fragmentation procedure. Im- 
mediately after the fragmentation the particles must overlap 
due to mismatch of geometry of the spherical fragments (mid- 
dle), a small repulsive force helps to overcome this unphysical 
transient situation. 



Let us remark that our algorithm was flexible enough 
to include multiple fracture events: By this we mean the 
observation that a fragment was apt to break itself even 
before its first complete separation from its original twin 
fragment. 

An animated sequ ence of fragmentation can be found 
in the Internet under http://summa.physik.hu-bcrlin.de/ 
~kies/mill/bm.html 
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Finally, we mention that the proposed algorithm is not 
the only way to model particle fragmentation. An alter- 
native algorithm has been put forward by Astrom and 
Herrmann EJ. 

3.2 Specification of the Model System 

As depicted in Fig. ||the cylindrical wall of our system has 
been modelled by an ensemble of spheres possessing the 
same material constants as the grains inside the container. 
However, the motion of wall spheres is not affected by 
impacts but, instead, strictly governed by the continuous 
rotation of the cylinder. Moreover, in spite of identical 
elastic and dissipative constants the wall particles were 
considered to be much harder so they never would break. 




Fig. 3. The cylindrical container is modeled by particles placed 
along the circumference of the container. Roughness is taken 
into account by stochastic variation of the sizes of these parti- 
cles. 

A suitable size distribution of wall particles generates 
the effect of surface roughness. This type of modelling 
has been used in several numerical simulations and is also 
used in experiments to guarantee well-defined boundary 
properties. The wall particles have been placed along the 
circumference of the cylinder. Just like in realistic ball 
mills we evenly distributed 16 toolbars along the inner 
periphery of the cylinder which serve to lift particles and 
thus prevent them from sliding downhill. We composed 
the toolbars again from a collection of rigid spheres. As 
reported in literature p2f the maximum milling efficiency 
is achieved for a filling degree, i.e. the ratio of material 
and container volume, of ca. 40%. 

The simulation results to be presented in the next sec- 
tion have been done with a system consiting in average 
of about 1000 particles. Due to fragmentation the num- 
ber of particles was not constant, i.e., the particle number 
fluctuates. The container was modeled by about 500 wall 
particles as described above. 

Initially the spheres which model the grist were ran- 
domly positioned inside the container avoiding any over- 
lap. This artificial initial condition has to be relaxed by 
evolution of the system. In this way we achieved a realistic 
configuration for any chosen rotation velocity. Obviously 



the typical asymptotic configuration depends on the ro- 
tation frequency v, therefore, we started our numerical 
monitoring only after transients had died out. In Fig. ^ 
we present snapshots of typical configurations depending 
on four different rotation frequencies. 

In technologically relevant ball mills there is a con- 
tinuous transport of raw material into the mill and of 
milled material out of the device. In most industrial ap- 
plications this exchange of raw material and final product 
is accomplished by axial transport (z-direction) of mate- 
rial. Of course, with our two-dimensional model we cannot 
account for axial transport, therefore, we simulated the 
material exchange in continuous operation mode in the 
following way: Whenever fragments were generated whose 
size dropped below a minimum we removed them from 
the container. After exact bookkeeping of the mass of re- 
moved "dust" particles we reinserted a big particle when 
the accumulated dust mass exceeded an upper threshold, 
of course, obeying mass conservation. In this way we con- 
served the average filling level throughout the simulation 
and, thus, could reach a steady state of the continuous 
operation mode. 



4 Simulation Results 

4.1 Mass distributions 

As a first result we studied the frequency dependent frag- 
ment size distribution. This distribution is not necessarily 
equal to the distribution of section which was the result 
of single fracture experiments. 

It is known pl| , p2p^ | that for an arbitrary fracture 
mechanism, which hierarchically is applied on all length 
scales, one will always generate a log-normal distribution 
of fragment sizes; this assertion is a pure consequence of 
the central limit theorem. However, in our modelling the 
ball mill is operated as an open system, i.e. there is a par- 
ticle flow into and out of the mill. Therefore, analyzing the 
asymptotic size distribution is also of theoretical interest. 

Furthermore, the explicit knowledge of the distribution 
plays also a practical role, because from that knowledge 
one can decide whether a multi-level fragmentation pro- 
cess is to be preferred to a single-level one or vice versa. 
A multi-level process can be either a distribution of the 
comminution process to several ball mills, where each is 
filled with granular material of a specified narrow size dis- 
tribution, or a spatial separation of material of different 
size in one mill. 

Fig. shows the differential fragment size distribution, 
i.e. the mass fraction of all particles with size within a cer- 
tain interval, for four different rotational frequencies. The 
black circles in each column correspond to measurements 
at different times and thus reflect the temporal variance. 
The time period of each simulation was 10 seconds real 
time (beyond transients). 
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Fig. 4. Typical dynamical configurations depend on the rotation speed: 0.5 Hz (top left), 1.0 Hz (top right), 1.25 Hz (bottom left), 
2.5 Hz. Grey scale codes the instantaneous maximum compression (normal force) experienced by each particle through contact 
with its neighbours (light: small normal force, dark: large force). Obviously the largest forces occur in the range of intermediate 
rotation velocities (1.0 Hz, 1.25 Hz). Moreover, we identify the location of largest forces which can initialize fragmentation to be 
close to the bottom of the rotating cylinder. This important observation will be discussed in detail in section H. 
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Fig. 5. The differential fragment size distribution for different 
rotation velocities. The black circles represent the values for a 
fixed mass interval at different times. The variance of the circles 
illustrates fluctuations of the distribution. The solid lines are 
tracing the averages of each interval. The left border of the 
scale is fixed by the size of the sieve (m s i e ve), the right border 
by the size of the refilled big grains (m raw ). 



Obviously the lowest curve (2.5 Hz) changes only very 
weakly during the simulation, which can be interpreted as 
a low rate of fragmentation. This curve - and with restric- 
tions also the top curve (0.5 Hz) - is determined by the 
initialization. Unfortunately, the middle curves are not sig- 
nificantly different, which mainly is a consequence of the 
short simulation time (10 seconds). To answer the inter- 
esting questions (theoretical distribution, multi-level pro- 
cess) it will be necessary to simulate considerably longer 
time periods. 

4.2 Optimization of the efficiency 

In the context of industrial applications the notion of ef- 
ficiency mainly focuses on the aspects of maximum com- 
minution rate and, perhaps, of power consumption. There 
are many free parameters which can be varied to maximize 
the efficiency: speed of rotation, filling degree, size distri- 
bution, cylinder diameter, grist size distribution, various 
mill types, etc. 

In our analysis we investigated the dependence of the 
comminution rate on the speed of rotation. In order to 
compare our numerical results with experimental data we 
normalized quantities in the following way: The speed of 
rotation is normalized to the velocity n c = yJg/M (2tt)~ 1 
for which the centrifugal force balances the gravitational 
force whereas the comminution rate is normalized to its 
measured maximum. The result of our simulation is plot- 
ted in Fig. H using filled circles and dashed line. Clearly, a 
maximum occurs around n « n c whereas the rate rapidly 
diminishes for larger or smaller velocities. 




0.5 1.0 1.5 

rel. speed of rotation n/n 



2.0 



Fig. 6. The dependence of the fraction rate on the rescaled 
rotation velocity, simulation (filled circles plus dashed line) and 
experiment (full line) taken from |}2| (original in |fjj| ). The 
error bars for the simulation visualize the fluctuations during 
the simulation period of approximately 10 seconds real time. 



The full line in Fig. ^ depicts a similar curve follow- 
ing from experiments. The curve shows the relative power 
consumption as a function of the normalized speed of ro- 
tation. Qualitatively both curves are in good agreement; 
systematic deviations can be explained by the following 
observation: Strictly speaking both curves relate to dif- 
ferent physical quantities. The fraction rate only accounts 
for the energy dissipated through fractures while the ex- 
perimental curve includes all dissipative processes, for in- 
stance friction between grains, frictions between the grains 
and the wall, dissipative impacts, heating, etc.). As a con- 
sequence, the curve for the simulation should be found 
below the experimental curve. Moreover, from the same 
reasoning it is also clear that observed deviations are mi- 
nor if the main portion of the overall dissipated energy 
results from the fractures, i.e. around the maximum. 



4.3 Preferred fragmentation locations 

The series of pictures in Fig. fj] illustrates the spatial distri- 
bution of fragmentation locations for four different rota- 
tion velocities. We restricted our presentation to the most 
important lower right segment of the ball mill (clockwise 
rotation) . To compute this distribution the coordinates of 
breaking particles were stored during the simulation. Af- 
ter the simulation the fraction of all fragmentation events 
occurring at places inside little cells (coarse graining) was 
computed. In Fig. [7] grey values code the frequency: dark 
means high frequency and light low frequency respectively. 

In the bottom right figure one finds a few fragmen- 
tation events near the cylinder wall. The restriction to 
medium- grey values results from the fact, that the rare 
events are equidistributed across the inner perimeter. 
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Fig. 7. Spatial distribution of fragmentation locations for different rotation velocities: 0.5 Hz, 1.0 Hz, 1.25 Hz, 2.5 Hz (from top 
left to bottom right). Grey values code the frequency of fragmentation in a certain region (cell), dark meaning high frequency 
and light low frequency respectively. Preferred fragmentation regions are found near the bottom of the mill, i.e. deep inside the 
granular material. This observation is in good agreement with experimental results and the results from chapter H. 
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The predominance of rather light grey (or even white) 
pixels in the rest of the figures hints at regions where 
fragmentation events are rare as compared to the sparse 
medium-grey to dark pixels which indicate regions of high 
fracture probability. The observation that most fragmen- 
tation events occur deep inside the material and not near 
the surface, which is heavily agitated by impacts, is, at 
first sight, somewhat counter-intuitive. An explanation for 
this fact will follow from a closer analysis of the local pres- 
sure concentration in the next section. 



5 Spatial distribution of pressure in ball mills 

In this chapter we will dwell on the physical origin of 
the unusual pressure distribution f45| measured in exper- 
iments. 

As can be concluded from Figs. fj] and [|, the major- 
ity of fragmentation events does not occur near the free 
surface of the granular material but deep inside the ma- 
terial near to the walls. This strange behavior was ob- 
served already in experiments by Rothkegel and Rolf p],|| 
in two-dimensional model mills. They used instrumented 
balls which emitted a flash of light whenever the strain ex- 
ceeded a fixed threshold. The positions of recorded flashes 
were digitized and from the statistics of these data they 
inferred the stationary strain distribution inside the ball 
mill. For various threshold values they always found the 
peak strain deep inside the material (see Fig. |J). This was 
surprising because inside the material relative grain veloc- 
ities are rather low in comparison with relative velocities 
at the surface. 

This experimental finding, together with its validation 
by our simulation (see the previous section), calls for an 
explanation. We will solve this challenging, so far open 
question using the powerful possibilities of molecular dy- 
namics, namely the option to monitor physical quantities 
which are hardly accessible through experiments. 

In the following simulations we ignored fragmentation 
events - which are not relevant for our explanation - sim- 
ply for the sake of computational efficiency, so we could 
afford longer simulation time. This time the radius of 
the two-dimensional ball mill was set to M = 4 cm. The 
mill was filled with N = 800 spherical grains with radii 
equally distributed in the interval [0, 05, 0, 11] cm. We per- 
formed simulations with three different rotation velocities 
J?i = 2 Hz, % = 10 Hz and 17 m = 19 Hz. 

Fig. ^ shows snapshots of our performed simulations. 
The pictures in the top row are taken form a simulation 
with rotation velocity 17 — Q\ . The velocity is sufficiently 
enough to observe a continuous flow at the surface. The 
surface shape is similar to a plane. The regime with dis- 
continuous flow (for lower rotation velocities) will not be 
investigated in this paper, because it is not relevant for 
technical purposes. For a detailed discussion of the in- 
teresting observations at the transition between continu- 
ous and stick-slip flow we refer to the work by Rajchen- 
bach ||. 

The picture in the middle row of Fig. || depicts snap- 
shots of simulation with 17 = i?n. The surface of the ma- 



terial no longer shapes a plain but the grains fall down a 
steep slope hitting a flat surface. 

The bottom row of Fig. ^ exhibits the snapshot of a 
simulation with the highest rotation velocity 17 = 17m. 
Here the grains are carried away with the fast rotation 
of the cylinder. Beyond a certain point they are likely 
to loose contact with the wall and follow a free parabola 
before impacting back down on the surface. 

Moreover, the snapshots of Fig. || provide an intu- 
itive understanding, what mechanism is responsible for 
the pressure distribution mentioned above. Again, grey 
values in these pictures codes the local pressure Pi acting 
on grain i given by 

P*=£if. ( 21 ) 
i 

The index j is running over all particles which are in con- 
tact with grain i. Dark circles denote high pressure, light 
spheres low pressure. 

Obviously many of the heavier stressed particles are 
deep inside the material, which is in good agreement with 
the experimental observation by Rothkegel and Rolf jl],|| . 
Most importantly, the pictures reveal that grains expe- 
riencing peak stress often form linear structures. In the 
following we focus on the development and the properties 
of such force chains. By numerical analysis we will prove 
that those force chains are the crucial physical reason for 
the observed pressure distribution. 

We defined force chains by a set of three self-evident 
conditions: Grains i, j and k are considered to be members 
of the same force chain if: 

1. Particles i, j and j, k are next neighbours. 

2. The pressure acting on each of the grains exceeds a 
certain threshold (Pi > 1000 gems -2 ). 

3. The connecting lines between i, j and j, k form an 
angle larger than 150°, i.e. the centres of three grains 
almost fall on a line. 

These three conditions are evaluated by a computer algo- 
rithm. In Fig. [)] all grains belonging to a force chain are 
graphically connected by lines. 

We see that most of the harder stressed grains could 
be assigned to a force chain. From this we conclude that 
the main part of the static and dynamic pressure prop- 
agates along the force chains. The pressure acting on a 
highly stressed grain inside a force chain is up to 100 times 
larger as compared to the average pressure of the neigh- 
bouring grains not belonging to a force chain, i.e. the force 
distribution is strongly inhomogeneous. 

The development of such force chains is not restricted 
to grains in a ball mill but has been observed in different 
granular systems by various experimentalists [ll].[T^.[lTi|.{^o| 
and numerically |5l|,^2|,^3) . It could be shown, that one of 
the granular phases is characterized by the observation 
of force chains [Q. From this one might conclude that 
the occurrence of force chains is an inherent feature of all 
granular matter, which may substantially contribute to its 
specific characteristics. 



10 



Volkhard Buchholtz, Jan A. Freund, Thorstcn Poschel: Molecular Dynamics of Comminution in Ball Mills 




Fig. 9. Snapshots of the simulation for different rotation velocities. Grey values code the local pressure Pi (see text). Linear 
segments connect grains belonging to a force chain. Qi — 2 Hz (top), Qu = 10 Hz (middle), J?m = 19 Hz (bottom). 
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is nothing but the number of particles belonging to the 
chain). These pairs are shown for three different rotation 
velocities varying from Q\ (top) to Q\u (bottom). 

The frequency of a force chain drops exponentially 
with increasing length L. For the highest rotation velocity 
(l?iu) the statistics becomes unreliable for chains longer 
than L = 13. This behaviour, which is in contrast to obser- 
vations for slower measured rotations, is due to the higher 
energy input, which causes a decompactification of the 
material. 

Concerning the earlier mentioned spatial strain distri- 
bution the maximal pressure measured in each force chain 
is of paramount interest rath er t han the average pressure 
in a force chain (see section 2.2). For lower rotation ve- 
locities f2\ and i?n a monotonous increase of the maximal 
pressure with the length L can be observed. For a dis- 
cussion of these curves one has to take into account that 
the pairs of plots in Fig. [l^ do not constitute indepen- 
dent quantities. Choosing a sample of L particles from an 
ensemble of particles on which forces with a given distri- 
bution are acting, the maximal force in this sample will in- 
crease with increasing L. Therefore, in order to prove sig- 
nificance, the result in Fig. has to be weighted against 
this purely statistical effect. 

The probability to measure a maximal pressure p max 
in the interval P max g [x; x + dx] in a sample of size L is 
equal to the probability that all particles feel a pressure 
P < x + dx minus the probability that all L particles feel 
a pressure P < x. 

p L (F ra e [x; x + dx}) = [p{P < x + dx)] L - [p{P < x)] L 

(22) 

Assuming equally distributed pressures P € [0 : P m ], 
p(P < x) means 



p(P < x) 



—d ' = — 
P P 



(23) 



By inserting Eq. ( |23| ) into Eq. 
order of dx 

Pi(r ax 6 [x;x + dx]) = 



one finds in first 



1 



(Pm) 1 



Lx L ~ 1 dx . 



(24) 



Fig. 10. Within each pair of the triplet (I,II,III) the lower 
plot presents the frequency distribution whereas the upper one 
shows the maximal pressure inside a force chain both as a 
function of its length L (abscissa for all plots). i?i = 2 Hz (top), 
Q\i = 10 Hz (center), Qui — 19 Hz (bottom) (explanation of 
lines see text). 



Fig. [To] comprises six plots which come in three pairs; 
the lower plot of each pair always depicts the frequency 
distribution of force chains whereas the related upper one 
presents the (statistically evaluated) maximal pressure both 
as a function of the length L of a force chain (which 



Therefore, the average maximal pressure as a function 
of sample size L is given by 



(P max (I)) = 



(Pm) 1 



xLx L 1 dx 



1 



L 



(PmT L + 1 

L 

. p 

1 m • 



-P. 



L+l 



L + l 



(25) 



The dashed line in Fig. [UJ shows this statistical effect for 
pmax w ^ p m — 4800gcms~ 2 . Obviously, these curves 
do not reach up to our numerical data of the top and 
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center figure; especially for larger L the dashed line falls 
far below. The numerical data points in the bottom figure 
(high rotation velocity) to some extent coincide with the 
dashed curve, which only represents the statistical effect 
discussed above. Hence, we can conclude that for lower 
velocities the maximum stress significantly concentrates 
within long force chains whereas for high rotation veloci- 
ties a significant effect of force chains is not evident. 

From experimental and numerical results (for instance 
p9| . p4|j5^ ]) it is known, that the observation of force chains 
is connected to a strongly inhomogeneous pressure dis- 
tribution. The distribution is essentially an exponentially 
decreasing function with increasing pressure. A simplified 
expression for this distribution is given by 



p{P) = a exp{— a P} 



(26) 



In analogy to Eq. ( p3| ) the probability to measure a 
pressure P < x is given by 

X 

p(P < x) = J a exp{— ax'}dx' = 1 — exp{-ax} . (27) 



Inserting this into Eq. (|22[) we obtain in the first order of 
dx 

PL (P max G [x;x + dx]) = 

I _ p -a(x+dx)^ L _ ^ _ g-aij 1 
L-l. 



L(l-e- ax ) L - L e- ax adx 



(28) 



The average of the maximal pressure can be calculated 
by integrating 



(F max (L)) = / xL(l-e- ax ) L 1 e- ax adx 



N 



= lim ix (\-e- ax ) ^ - / (1 -e~ ax ) dx } (29) 



An iterative splitting 

(i - e - ax y = (i - e - ax y- 1 (i - e~ ax ) 

results in a summation of integrals of the type 

N 

l^e- ax ) l - 1 e- ax dx^—. 

ai 

o 



Thus, we obtain 



(P ma %L))=]T- 
L — ' at 

i—l 



(30) 



(31) 



(32) 



A restriction of the statistics to values of the pressure 
P > 1000 g cm s -2 , performing an analogous calculation, 
results in 



L 1 

(P max (L)) = V — + 1000 gems 
t— i a i 



From the simulations we extracted the pressure distribu- 
tion yielding - as expected - an exponential decrease. The 
prefactor is a = 0.0009 for Q = 10 Hz and a = 0.001 for 
[2 = 2 Hz. The corresponding curves are plotted in Fig. [Ic] 
as dotted lines and reasonably describe to the measured 
values. Thus, we indirectly can conclude from the mea- 
surement that we have an exponential distribution, which 
indicates the occurrence of force chains. 
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Fig. 11. The maximal force acting on a particle as a function 
of the relative position in the force chain. The monotonous be- 
havior indicates that forces are accumulating downwards along 
the chain. 



Further insight in the mechanism of force chains is 
achieved through Fig. [n] which shows the average of the 
maximal force acting on a grain as a function of its posi- 
tion within the force chain. The average was performed for 
all force chains of length L = 8. The top particle of a chain, 
corresponding to the maximal y-coordinate, was defined 
as position 1. The monotonous increase of this maximal 
force indicates that, on the average, strain accumulates 
downwards along the force chain. In this respect a force 
chain is similar to a beam in a framework: The masses and 
momenta of all particles above a force chain are supported 
by this chain, with the consequence that lower particles of 
the chain are heavier loaded than upper ones. This trend is 
clearly visible in the data of Fig. O. Furthermore, a force 
chain can shield neighbouring particles, i.e. these grains 
experience a much weaker force which can be observed in 
the simulation data as well. 

Figure [l^ shows the spatial pressure distribution for 
different rotation velocities. Here grey values code the ab- 
solute local pressure, with light pixels indicating low and 
dark pixels high values, in order to allow for a direct com- 
parison of all three images. 

The presented data are temporal averages. To extract 
the spatial distribution of this data, the two-dimensional 
space was coarse-grained. For each time step the positions 
of all particles were mapped to the grid and the result- 
ing instantaneous local pressure was added to its related 
cell. Notice that the averaged data does not reflect the 
occurrence of force chains. 
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Fig. 12. Spatial distribution of the pressure. Grey values code the absolute pressure, with light pixels indicating low and dark 
pixels high pressure values. The grey scale is unique for all three figures to allow for a direct comparison. From left to right: 
Q — 2 Hz, Q = 10 Hz, O = 19 Hz. While for lower rotation velocities the maximum is deep inside the material near to the wall, 
for higher rotation velocity it is near to the free surface, which is heavily agitated by impacts. 



> ■ 




Fig. 13. Spatial frequency distribution of particles, exerted to a pressure P > 3000 gems -2 . Grey values, unique throughout 
the images, code the frequency: dark pixels high frequency and light pixels low frequencies. The top row only accounts for 
particles which do not belong to a force chain while, in contrast, the bottom row includes only particles being member of a force 
chain. From left to right: O = 2 Hz, O = 10 Hz, O = 19 Hz. By visual inspection it becomes obvious that most of the heavier 
loaded particles belong to force chains and are found near to the wall of the mill. 
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Surprisingly, for lower rotation velocities one finds that 
direct impacts with high relative velocities at the material 
surface do not result in high local pressure. Therefore, 
these impacts are not relevant for the fragmentation of 
grains in ball mills. Regions of high pressure can be found 
mainly near to the wall, deep inside the material. This 
observation is in good agreement with the experimental 
result of Rothkegel and Rolf [QJD . The extreme right figure 
in Fig. [[2] reveals a different behavior. The maximum is 
near to the material surface where heavy impacts occur. 
The absolute values of the pressure are much lower than 
in the left figures. This observation relates to the existence 
of rather weak force chains. 

Figs. |l3| show the spatial distribution of the number of 
particles, which experience a pressure Pj > 3000 g cm s -2 . 
In these figures the particles which belong to a force chain 
are only considered in the bottom row whereas particles 
not belon ging to a force chain are contained in the top 
row. Figs. 13| can be directly compared with the figures of 
Rothkegel |M (see Fig. ||); both are found to be in good 
agreement. The maxima of the distribution for fl\ and 
flu lie near to the wall of the mill, inside the material 
as observed in the experiment. Those maxima are an ex- 
clusive result of force chains. Grains not belonging to a 
force chain contribute only weakly to the comminution. 
For the higher velocity Dm direct impacts of particles at 
the surface are the dominant part. Nevertheless, absolute 
values are much lower than for i?i and Q\\. Again, this 
observation is in good agreement with the experiment jj| . 

6 Conclusion 

We have simulated the process of autogenous dry com- 
minution in a ball mill by the method of molecular dy- 
namics. To this end we developed a novel algorithm which 
accounts for fragmentation of particles. Throughout our 
simulations all particles, including the fragments, were 
modelled as spheres. To accomplish this idealization we 
introduced rather short-lived unphysical transient situa- 
tions in which the fragments resulting from a cracked par- 
ticle penetrate each other. The results of our simulations 
are available in the internet as animated sequences und er 
http: / /summa.physik.hu-bcrlin.dc/^kies/mill/bm.html 

Employing our molecular dynamics algorithm we could 
satisfactorily reproduce experimental phenomena, e.g. the 
normalized comminution rate as a function of the normal- 
ized rotation speed and the spatial pressure distribution. 

From the achieved numerical data we were able to ex- 
plain an experimental result which was poorly understood 
so far : Comminution processes in ball mills occur deep 
inside the material rather than close to the surface where 
the relative collision velocity is maximal. We explained 
this effect by the formation of force chains. Particles which 
are members of one ore more force chains experience a sig- 
nificantly larger average force, thus substantially enhanc- 
ing the fragmentation probability, as compared to parti- 
cles which do not belong to any force chain. 

From this we concluded that the experimentally ob- 
served and numerically achieved spatial pressure and frag- 



mentation distributions are crucially determined by the 
formation of force chains. Direct collisions of particles are 
of minor importance for comminution. Without the ex- 
istence of force chains ball mills would work much less 
efficient. In the context of a ball mill efficiency optimiza- 
tion it would be intriguing to investigate the possibility 
of enhancing the formation of force chains, e.g. by modi- 
fications of the geometrical properties of mills, additional 
application of vibration etc. 

The method developed in this article may be elabo- 
rated in many different aspects: 

— The numerical investigations were done in two dimen- 
sions. Ball mills are three-dimensional devices. 

— Axial transport of the material was not simulated. 

— All particles, including fragments, were spherical. 

— Only a single-level process was simulated. Modern ball 
mills frequently operate in multi-level mode. 

— Adhesive forces were not considered in the force law. 

— Variation of material properties which originates from 
an increase of temperature due to comminution was 
not considered. 

We expect that molecular dynamics will rapidly de- 
velop towards an important tool for the powder technol- 
ogy machinery. Due to the rapid progress in hardware and 
software technology, including parallel algorithms, it will 
be possible soon to construct and optimize powder tech- 
nology facilities. 
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